# remove.packages("brainMapR")
devtools::install_github("baptisteCD/brainMapR")
library(brainMapR)

system.file("extdata/07_BWAS_MOA_realPheno", "", package = "brainMapR",mustWork = TRUE)

system.file("extdata/03_BWAS_uncorrected_realPheno", "", package = "brainMapR",mustWork = TRUE)

for (pheno in c("Age_MRI","body_mass_index_bmi_f21001_2_0","fluid_intelligence_score_f20016_2_0", "sexuses_datacoding_9_f31_0_0","smoking_status_f20116_2_0", "sexuses_datacoding_9_f31_0_0") ){
for (scenario in c("03_BWAS_uncorrected_realPheno", "09_BWAS_ICVagesex_realPheno", "04_BWAS_5globalPCs_realPheno", "05_BWAS_10globalPCs_realPheno", "06_BWAS_10specificPCs_realPheno", "07_BWAS_MOA_realPheno",  "10_BWAS_MOA_multiORM_QC_realPheno", "15_BWAS_MOA_allModa_FE_realPheno")){
  
identifyClustersBWAS(pathToTheBWASresults =system.file(paste0("extdata/", scenario), "", package = "brainMapR",mustWork = TRUE) , bwasFile = paste0("BWAS_",pheno), outputFolder = "../", signifThreshold = 5e-8)
  print(paste(pheno, scenario))
}
}

cols=c(RColorBrewer::brewer.pal(n = 10, name = "RdYlBu")[10:6],RColorBrewer::brewer.pal(n = 10, name = "RdYlBu")[5:1]) # Select palette colours

png("legendbar.png", width = 5, height = 15, units = "cm", res = 400)
par(mar=c(2,4,2,2))
my.colors = colorRampPalette(cols)
z=matrix(1:100,nrow=1)
x=1
y=seq(-0.4,0.4,len=100) # supposing 3 and 2345 are the range of your data
image(x,y,z,col=my.colors(20),axes=FALSE,xlab="",ylab="", main="Correlation")
axis(2)
dev.off()

library(knitr)
include_graphics(path = "examplePlots/legendbar.png", dpi = 50)

library(viridis)
cols=viridis_pal(option = "C")(7)

for (var in c("Age_MRI","body_mass_index_bmi_f21001_2_0","fluid_intelligence_score_f20016_2_0", "sexuses_datacoding_9_f31_0_0","smoking_status_f20116_2_0", "sexuses_datacoding_9_f31_0_0") ){
  for (scenario in c("03_BWAS_uncorrected_realPheno", "09_BWAS_ICVagesex_realPheno", "04_BWAS_5globalPCs_realPheno", "05_BWAS_10globalPCs_realPheno", "06_BWAS_10specificPCs_realPheno", "07_BWAS_MOA_realPheno",  "10_BWAS_MOA_multiORM_QC_realPheno", "15_BWAS_MOA_allModa_FE_realPheno")){

plotCortical(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../")

plotSubcortical(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../")

plotSubcortical_flat(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../")

combineCorticalSubcorticalPlots(inputPath = "../", bwasFile = paste0("BWAS_",pheno), pathPhenotypeFile = "../../../UKB_phenotypes15K/Age_MRI.phen", outputPath = " ../" )
}
}

library(viridis)
cols=viridis_pal(option = "C")(7)

BrainMapManhattanPlot(inputPath = system.file("extdata/07_BWAS_MOA_realPheno", "", package = "brainMapR",mustWork = TRUE), path = "03_BWAS_uncorrected_realPheno/", bwasFile = "BWAS_Age_MRI.linear", yMax = 120, signifThreshold =5e-8,  phenotypeLabel = "Age")

include_graphics(path = paste0("examplePlots/03_BWAS_uncorrected_realPheno/Manhathan_BWAS_Age_MRI.linear_simple.png"))

pheno=read.table("UKB_BWASvars_labels.txt", header=T, stringsAsFactors = F)
pheno=pheno[which(pheno$variable %in% c( "body_mass_index_bmi_f21001_2_0","fluid_intelligence_score_f20016_2_0", "Age_MRI","sexuses_datacoding_9_f31_0_0", "smoking_status_f20116_2_0")),]

vall=NULL
for (iii in 1:5){
  all=NULL
  for (scenario in c("03_BWAS_uncorrected_realPheno","09_BWAS_ICVagesex_realPheno", "05_BWAS_10globalPCs_realPheno", "07_BWAS_MOA_realPheno", "15_BWAS_MOA_allModa_FE_realPheno", "10_BWAS_MOA_multiORM_QC_realPheno")){
  
  vres=read.csv(paste0("../../",scenario, "/", "BWAS_signif_", pheno$variable_clean[iii], ".csv" ))
  res=read.table(paste0("../../",scenario, "/", "Results_clusterFWER_", pheno$variable[iii], ".txt" ), header=T)
  nclust=res[,c("NumberClusters_thickness", "NumberClusters_area", "NumberClusters_thick", "NumberClusters_LogJacs")]
  maxclust=res[,c("maxSizeCluster_thickness", "maxSizeCluster_area", "maxSizeCluster_thick", "maxSizeCluster_LogJacs")]
  
  all=rbind(all,  dim(vres)[1], sum(nclust), max(maxclust, na.rm = T))

  }
  vall=cbind(vall, all)
}

colnames(vall)=pheno$variable_clean
writexl::write_xlsx(as.data.frame(vall), "../../Summary_BWAS_realPheno_bonferroniNtraits_R1IEEE.xlsx")


# N clusters cortex
vall=NULL
for (iii in c(1,2,3,5)){
  all=NULL
  for (scenario in c("03_BWAS_uncorrected_realPheno","09_BWAS_ICVagesex_realPheno", "05_BWAS_10globalPCs_realPheno", "07_BWAS_MOA_realPheno","15_BWAS_MOA_allModa_FE_realPheno", "10_BWAS_MOA_multiORM_QC_realPheno")){
  #vres=read.csv(paste0("../../",scenario, "/", "BWAS_signif_", pheno$variable_clean[iii], ".csv" ))
  res=read.table(paste0("../../",scenario, "/", "Results_clusterFWER_", pheno$variable[iii], ".txt" ), header=T)
  nclust=res[,c("NumberClusters_thickness", "NumberClusters_area", "NumberClusters_thick", "NumberClusters_LogJacs")]
  maxclust=res[,c("maxSizeCluster_thickness", "maxSizeCluster_area", "maxSizeCluster_thick", "maxSizeCluster_LogJacs")]
  all=rbind(all, sum(nclust[,1:2]))
  }
  vall=cbind(vall, all)
}
colnames(vall)=pheno$variable_clean

# Brain plot 
source("RR_9.1_BWAS_manhathanPlot_otherFunctions.R")
library(Morpho)
library(plyr)
library(viridis)


for (scenario in c( "07_BWAS_MOA_realPheno")){
# Plot cortical and subcortical
plotCortical_noSignif(path = paste0("path/to/working/directory/", scenario, "/"), variable = "body_mass_index_bmi_f21001_2_0")
plotSubcortical_noSignif(path =paste0("path/to/working/directory/", scenario, "/"), variable = "body_mass_index_bmi_f21001_2_0")
combineCorticalSubcorticalPlots_noSignif(path = paste0("path/to/working/directory/", scenario, "/"), variable = "body_mass_index_bmi_f21001_2_0")
}